** Regression analysis of EI across 4 waves
** Uses data generated from eisurveys_dataorg.do and then modified in fourwaves_descriptives to create additional measures
clear
set scheme s1color 

** Open monthly data
use thirteenmonths_analysis, clear

** Code model variables
gen under5=0 if hh_under5==0
replace under5=1 if hh_under5>=1 & hh_under5!=.

gen disab=0 if hh_disability==0
replace disab=1 if hh_disability>=1 & hh_disability!=.

gen over_150_fpl=0 if fpl!=.
replace over_150_fpl=1 if fpl==3

gen between_100_150_fpl=0 if fpl!=.
replace between_100_150_fpl=1 if fpl==2

gen lower_100_fpl=0 if fpl!=.
replace lower_100_fpl=1 if fpl==1

gen other=1 if race>3 & race!=1
replace other=0 if other==.

replace dwelling_condition=1 if dwelling_condition>0
rename dwelling_condition housing_conditions

gen asf_home=1 if dwelling_type==4
replace asf_home=0 if dwelling_type!=4 & dwelling_type!=.

gen dsf_home=1 if dwelling_type==3
replace dsf_home=0 if dwelling_type!=3 & dwelling_type!=.

gen mobile=1 if dwelling_type==1
replace mobile=0 if dwelling_type!=1 & dwelling_type!=.

gen apartment=1 if dwelling_type==2
replace apartment=0 if dwelling_type!=2 & dwelling_type!=.

gen otherhome=1 if dwelling_type==5
replace otherhome=0 if dwelling_type!=5 & dwelling_type!=.

** Combine other dwelling with apartment or condominium
replace apartment=1 if otherhome==1

gen payment_rent=1 if fuel_payment==2 & month<5
replace payment_rent=1 if fuel_payment==3 & month<5
sort caseid_w1 month
replace payment_rent=1 if fuel_payment[_n-1]==2 & caseid_w1==caseid_w1[_n-1] & month>4
replace payment_rent=1 if fuel_payment[_n-1]==3 & caseid_w1==caseid_w1[_n-1] & month>4
replace payment_rent=1 if fuel_payment[_n-2]==2 & caseid_w1==caseid_w1[_n-2] & month>4
replace payment_rent=1 if fuel_payment[_n-2]==3 & caseid_w1==caseid_w1[_n-2] & month>4
replace payment_rent=0 if payment_rent==.

** Generate alternative measure of FPL, that includes respondents with changes in income to above 200%+ FPL
gen fpl2=fpl
replace fpl2=4 if fpl==.
label define fpl 1 "100% FPL or below" 2 "100%-150% FPL" 3 "150%-200% FPL" 4 "200% FPL or above"
label values fpl2 fpl
gen under100fpl=1 if fpl2==1
replace under100fpl=0 if fpl2!=1
gen between100150fpl=1 if fpl2==2
replace between100150fpl=0 if fpl2!=2
gen between150200fpl=1 if fpl2==3
replace between150200fpl=0 if fpl2!=3
gen over200fpl=1 if fpl2==4
replace over200fpl=0 if fpl2!=4

** Move people with FPL over 200 into FPL 150-200
** These are people whose income increased from time of W1 screening to later in the survey
replace between150200fpl=1 if over200fpl==1

xi i.state
xi i.wave
xi i.month

recode hhbill notice disconnect(100=1)

*global demographics "age female black hispanic other under100fpl between100150fpl educ hh_under5 hh_disability electronic_device retired unemployed housing_conditions payment_rent pcturban"

global demographics "age female black hispanic other under100fpl between100150fpl educ hh_under5 hh_disability electronic_device retired unemployed housing_conditions"

global dwelling "dsf_home mobile apartment"

global govassist "wap liheap gov_assistance_year"

** Add labels for tables/figures
lab var hhbill "Unable to pay bill"
lab var notice "Received shut-off notice"
lab var disconnect "Disconnected from service"
lab var hh_under5 "Children age 5 years and younger"
lab var hh_disability "Disabled person in household"
lab var unemployed "Unemployed"
lab var electronic_device "Electronic medical device"
lab var retired "Retired"
lab var black "Black"
lab var hispanic "Hispanic" 
lab var other "Other nonwhite race/ethnicity"
lab var female "Female"
lab var age "Age"
lab var under100fpl "Household income below 100% FPL"
lab var between100150fpl "Household income between 100-150% FPL"
lab var between150200fpl "Household income between 150-200% FPL"
lab var educ "Educational attainment"
lab var asf_home "Attached single family home"
lab var dsf_home "Detached single family home"
lab var mobile "Mobile home or trailer"
lab var apartment "Apartment, condominium, or other"
lab var otherhome "Other dwelling type"
lab var housing_conditions "Housing conditions"
lab var gov_assistance_year "Government assistance last year"
lab var liheap "Received LIHEAP last month"
lab var wap "Received WAP last month"

** Table S1 (more below for additional DVs)
svy: mean hhbill notice disconnect age female black hispanic other white under100fpl between100150fpl between150200fpl educ hh_under5 hh_disability electronic_device retired unemployed housing_conditions asf_home dsf_home mobile apartment wap liheap gov_assistance_year

** Bill Pay
logit hhbill $demographics $dwelling $govassist i.state i.month [pweight=weight], vce(cl state)
estimates store bill 
** Table S2, column 1
outreg2 using table_s2, cti(Unable to pay bill) replace bdec(2) 10pct drop(i.state) addstat(Pseudo R-squared, `e(r2_p)') eform

** Predicted probabilities
** https://stats.oarc.ucla.edu/stata/dae/using-margins-for-predicted-probabilities/
logit hhbill age female i.black i.hispanic other i.under100fpl between100150fpl educ i.hh_under5 hh_disability i.electronic_device retired unemployed i.housing_conditions dsf_home mobile apartment i.gov_assistance_year i.liheap i.wap i.state i.month [pweight=weight], vce(cl state)
margins black hispanic under100fpl hh_under5 electronic_device housing_conditions gov_assistance_year liheap wap, atmeans

** Table S4
** Month-by-month
logit hhbill $demographics $dwelling $govassist [pweight=weight] if month==1, vce(cl state)
outreg2 using table_s4, cti(Month 1) replace bdec(2) 10pct drop(i.state) addstat(Pseudo R-squared, `e(r2_p)') eform label

forvalues i = 2(1)13 {
	logit hhbill $demographics $dwelling $govassist [pweight=weight] if month==`i', vce(cl state)
	estimate store month`i'
	outreg2 using table_s4, cti(Month`i') append bdec(2) 10pct drop(i.state) addstat(Pseudo R-squared, `e(r2_p)') eform label
}

** Notice
logit notice $demographics $dwelling $govassist i.state i.month [pweight=weight], vce(cl state) 
estimates store notice
** Table S2, column 2
outreg2 using table_s2, cti(Received Shut-off Notice) append bdec(2) 10pct drop(i.state) addstat(Pseudo R-squared, `e(r2_p)') eform label
estimates store notice

** Table S5
** Month-by-month
logit notice $demographics $dwelling $govassist [pweight=weight] if month==1, vce(cl state) 
outreg2 using table_s5, cti(month`i') replace bdec(2) 10pct drop(i.state) addstat(Pseudo R-squared, `e(r2_p)') eform label

forvalues i = 2(1)13 {
	logit notice $demographics $dwelling $govassist [pweight=weight] if month==`i', vce(cl state) 
	outreg2 using table_s5, cti(month`i') append bdec(2) 10pct drop(i.state) addstat(Pseudo R-squared, `e(r2_p)') eform label
}

** Disconnect
** Pooling across months
logit disconnect $demographics $dwelling $govassist i.state i.month [pweight=weight], vce(cl state)
estimates store disconnect 
** Table S2, column 3
outreg2 using table_s2, cti(Disconnected from Service) append bdec(2) 10pct addstat(Pseudo R-squared, `e(r2_p)') eform label

** Table S6
** Month-by-month
logit disconnect $demographics $dwelling $govassist [pweight=weight] if month==1, vce(cl state) 
outreg2 using table_s6, cti(month`i') replace bdec(2) 10pct addstat(Pseudo R-squared, `e(r2_p)') eform label

forvalues i = 2(1)13 {
	logit disconnect $demographics $dwelling $govassist [pweight=weight] if month==`i', vce(cl state) 
	outreg2 using table_s6, cti(month`i') append bdec(2) 10pct addstat(Pseudo R-squared, `e(r2_p)') eform label
}

** Figure 3
** Coefplot of estimates from bill, notice, and disconnect
coefplot (bill, label(Unable to pay bill)) (notice, label(Received shut-off notice) msymbol(T)) (disconnect, label(Disconnected from service) msymbol(S)), keep(age female black hispanic other under100fpl between100150fpl educ hh_under5 electronic_device retired unemployed housing_conditions dsf_home mobile apartment otherhome gov_assistance_year liheap wap) xline(1) legend(symxsize(2) row(1)) msize(tiny) ylabel(,labsize(vsmall)) xlabel(,labsize(vsmall)) xlabel(-2(1)25) legend(size(vsmall)) eform saving(figure3, replace)


** Table S3
** OLS model
reg hhbill $demographics $dwelling $govassist i.state i.month [pweight=weight], cluster(state) 
outreg2 using table_s3, cti(Unable to pay bill) replace bdec(2) 10pct drop(i.state) label
reg notice $demographics $dwelling $govassist i.state i.month [pweight=weight], cluster(state) 
outreg2 using table_s3, cti(Received shut-off notice) append bdec(2) 10pct drop(i.state) label
reg disconnect $demographics $dwelling $govassist i.state i.month [pweight=weight], cluster(state)
outreg2 using table_, cti(Disconnected from service) append bdec(2) 10pct drop(i.state) label


** Analysis of progression/persistence

** Generate lagged DVs for one, two, and three months
sort caseid_w1 month
gen hhbill_lag1 = hhbill[_n-1] if caseid_w1==caseid_w1[_n-1]
gen hhbill_lag2 = hhbill[_n-2] if caseid_w1==caseid_w1[_n-2]
gen hhbill_lag3 = hhbill[_n-3] if caseid_w1==caseid_w1[_n-3]

gen notice_lag1 = notice[_n-1] if caseid_w1==caseid_w1[_n-1]
gen notice_lag2 = notice[_n-2] if caseid_w1==caseid_w1[_n-2]
gen notice_lag3 = notice[_n-3] if caseid_w1==caseid_w1[_n-3]

gen disconnect_lag1 = disconnect[_n-1] if caseid_w1==caseid_w1[_n-1]
gen disconnect_lag2 = disconnect[_n-2] if caseid_w1==caseid_w1[_n-2]
gen disconnect_lag3 = disconnect[_n-3] if caseid_w1==caseid_w1[_n-3]

lab var hhbill_lag1 "Unable to pay bill last month"
lab var notice_lag1 "Received shut-off notice last month"
lab var disconnect_lag1 "Disconnected from service last month"

** Generate lagged, 3-month dichotomous measure
sort caseid_w1 month
gen hhbill_lag=1 if hhbill_lag1==1 | hhbill_lag2==1 | hhbill_lag3==1
replace hhbill_lag=0 if hhbill_lag1==0 & hhbill_lag2==0 & hhbill_lag3==0

gen notice_lag=1 if notice_lag1==1 | notice_lag2==1 | notice_lag3==1
replace notice_lag=0 if notice_lag1==0 & notice_lag2==0 & notice_lag3==0

gen disconnect_lag=1 if disconnect_lag1==1 | disconnect_lag2==1 | disconnect_lag3==1
replace disconnect_lag=0 if disconnect_lag1==0 & disconnect_lag2==0 & disconnect_lag3==0

** Figure 5 and Table S10
** Create figure with just lagged terms for all three EI indicators
logit hhbill $demographics $dwelling $govassist hhbill_lag1 notice_lag1 disconnect_lag1 i.state i.month [pweight=weight], vce(cl state)
outreg2 using table_s10, cti(Bill) replace bdec(2) 10pct keep($demographics $dwelling $govassist hhbill_lag1 notice_lag1 disconnect_lag1) label eform
estimates store bill

logit notice $demographics $dwelling $govassist hhbill_lag1 notice_lag1 disconnect_lag1 i.state i.month [pweight=weight], vce(cl state)
outreg2 using table_s10, cti(Notice) append bdec(2) 10pct keep($demographics $dwelling $govassist hhbill_lag1 notice_lag1 disconnect_lag1) label eform
estimates store notice

logit disconnect $demographics $dwelling $govassist hhbill_lag1 notice_lag1 disconnect_lag1 i.state i.month [pweight=weight], vce(cl state)
outreg2 using table_s10, cti(Disconnect) append bdec(2) 10pct keep($demographics $dwelling $govassist hhbill_lag1 notice_lag1 disconnect_lag1) label eform
estimates store disconnect

** Coefplot of estimates from bill, notice, and disconnect
coefplot (bill, label(Unable to pay bill)) (notice, label(Received shut-off notice) msymbol(T)) (disconnect, label(Disconnected from service) msymbol(S)), keep(hhbill_lag1 notice_lag1 disconnect_lag1) xline(1) legend(symxsize(2) row(1)) msize(small) ylabel(,labsize(vsmall)) xlabel(,labsize(vsmall)) xlabel(-2(1)20) legend(size(vsmall)) eform saving(figure5, replace)
estimates clear

** Summary stats of lagged measures for Table 2
svy: mean hhbill_lag1 notice_lag1 disconnect_lag1

** Figure S6
** Robustness -- any time in the last three months
** Create figure with just lagged terms for all three EI indicators
logit hhbill $demographics $dwelling $govassist hhbill_lag notice_lag disconnect_lag i.state i.month [pweight=weight], vce(cl state)
estimates store bill

logit notice $demographics $dwelling $govassist hhbill_lag notice_lag disconnect_lag i.state i.month [pweight=weight], vce(cl state)
estimates store notice

logit disconnect $demographics $dwelling $govassist hhbill_lag notice_lag disconnect_lag i.state i.month [pweight=weight], vce(cl state)
estimates store disconnect

** Coefplot of estimates from bill, notice, and disconnect
lab var hhbill_lag "Unable to pay bill during last 3 months"
lab var notice_lag "Received shut-off notice during last 3 months"
lab var disconnect_lag "Disconnected from service during last 3 months"

coefplot (bill, label(Unable to pay bill)) (notice, label(Received shut-off notice) msymbol(T)) (disconnect, label(Disconnected from service) msymbol(S)), keep(hhbill_lag notice_lag disconnect_lag) xline(1) legend(symxsize(2) row(1)) msize(small) ylabel(,labsize(vsmall)) xlabel(,labsize(vsmall)) xlabel(-2(2)40) legend(size(vsmall)) eform 


** Analysis of Short-Term (Acute) v. Persistent EI (Chronic)
estimates clear

** Create categories of EI for ordered logit models
gen hhbill_cat=0 
replace hhbill_cat=1 if hhbill_acute==1
replace hhbill_cat=2 if hhbill_chronic==1
label define bill_chronic 0 "No Bill EI" 1 "Short-term" 2 "Persistent" 
label values hhbill_cat bill_chronic

gen notice_cat=0 
replace notice_cat=1 if notice_acute==1
replace notice_cat=2 if notice_chronic==1
label define notice_chronic 0 "No Notice EI" 1 "Short-term" 2 "Persistent" 
label values notice_cat notice_chronic

gen disconnect_cat=0 
replace disconnect_cat=1 if disconnect_acute==1
replace disconnect_cat=2 if disconnect_chronic==1
label define disconnect_chronic 0 "No Disconnect EI" 1 "Short-term" 2 "Persistent" 
label values disconnect_cat disconnect_chronic

** Create categories for successive month analysis
gen bill_atleast2=1 if hhbill_2months==1 | hhbill_3months==1 | hhbill_4months==1 | hhbill_5months==1 | hhbill_6months==1
replace bill_atleast2=0 if bill_atleast2==. 

gen notice_atleast2=1 if notice_2months==1 | notice_3months==1 | notice_4months==1 | notice_5months==1 | notice_6months==1
replace notice_atleast2=0 if notice_atleast2==. 

gen disconnect_atleast2=1 if disconnect_2months==1 | disconnect_3months==1 | disconnect_4months==1 | disconnect_5months==1 | disconnect_6months==1
replace disconnect_atleast2=0 if disconnect_atleast2==. 

** Only estimate regressions for respondents that participated in all four waves
** Collapse to Wave 4 respondents only
sort caseid_w1
by caseid_w1: egen wave_max=max(wave)
keep if wave_max==4

collapse (max) hhbill* bill* notice* disconnect* age female black white hispanic other under100fpl between100150fpl between150200fpl educ hh_under5 hh_disability electronic_device retired unemployed housing_conditions asf_home dsf_home mobile apartment wap liheap gov_assistance_year state wave month weight, by(caseid_w1)

xi i.state

** Add labels for tables/figures
lab var hhbill "Unable to pay bill"
lab var notice "Received shut-off notice"
lab var disconnect "Disconnected from service"
lab var hh_under5 "Children age 5 years and younger"
lab var hh_disability "Disabled person in household"
lab var unemployed "Unemployed"
lab var electronic_device "Electronic medical device"
lab var retired "Retired"
lab var black "Black"
lab var hispanic "Hispanic" 
lab var other "Other nonwhite race/ethnicity"
lab var female "Female"
lab var age "Age"
lab var under100fpl "Household income below 100% FPL"
lab var between100150fpl "Household income between 100-150% FPL"
lab var between150200fpl "Household income between 150-200% FPL"
lab var educ "Educational attainment"
lab var asf_home "Attached single family home"
lab var dsf_home "Detached single family home"
lab var mobile "Mobile home or trailer"
lab var apartment "Apartment, condominium, or other"
lab var housing_conditions "Housing conditions"
*lab var pcturban "Urban Percent"
*lab var payment_rent "Pays energy bill in rent"
lab var gov_assistance_year "Government assistance last year"
lab var liheap "Received LIHEAP last month"
lab var wap "Received WAP last month"

** Table 2 (summary statistics of additional DVs)
svy: mean hhbill_cat notice_cat disconnect_cat hhbill_number notice_number disconnect_number bill_atleast2 notice_atleast2 disconnect_atleast2


** Analysis of ordered categories of EI (ordered logit models)
** Bill Pay
ologit hhbill_cat $demographics $dwelling $govassist i.state [pweight=weight], vce(cl state) or
estimates store bill
** Table S7, Column 1
outreg2 using table_s7, cti(Unable to pay bill) replace bdec(2) 10pct drop(i.state) addstat(Pseudo R-squared, `e(r2_p)') eform label

** Predicted probabilities
ologit hhbill_cat age i.female i.black i.hispanic i.other i.under100fpl i.between100150fpl i.educ i.hh_under5 i.hh_disability i.electronic_device i.retired i.unemployed i.housing_conditions i.dsf_home i.mobile i.apartment i.wap i.liheap i.gov_assistance_year i.state [pweight=weight], vce(cl state) or

margins, dydx(black hispanic hh_under5 electronic_device housing_conditions wap) predict(outcome(0)) 
** Table 1, Column 1
margins, dydx(black hispanic hh_under5 electronic_device housing_conditions wap) predict(outcome(1))
** Table 1, Column 2
margins, dydx(black hispanic hh_under5 electronic_device housing_conditions wap) predict(outcome(2))

** Test of parallel lines assumption
omodel logit hhbill_cat age female black hispanic other under100fpl between100150fpl educ hh_under5 hh_disability electronic_device retired unemployed housing_conditions dsf_home mobile apartment wap liheap gov_assistance_year

** Notice
ologit notice_cat $demographics $dwelling $govassist i.state [pweight=weight], vce(cl state) or
estimates store notice
** Table S7, Column 2
outreg2 using table_s7, cti(Received a shut-off notice) append bdec(2) 10pct drop(i.state) addstat(Pseudo R-squared, `e(r2_p)') eform label

** Predicted probabilities
ologit notice_cat age i.female i.black i.hispanic i.other i.under100fpl i.between100150fpl i.educ i.hh_under5 i.hh_disability i.electronic_device i.retired i.unemployed i.housing_conditions i.dsf_home i.mobile i.apartment i.wap i.liheap i.gov_assistance_year i.state [pweight=weight], vce(cl state) or

margins, dydx(black hispanic hh_under5 electronic_device housing_conditions wap) predict(outcome(0)) 
** Table 1, Column 3
margins, dydx(black hispanic hh_under5 electronic_device housing_conditions wap) predict(outcome(1))
** Table 1, Column 4
margins, dydx(black hispanic hh_under5 electronic_device housing_conditions wap) predict(outcome(2))

** Test of parallel lines assumption
omodel logit notice_cat age female black hispanic other under100fpl between100150fpl educ hh_under5 hh_disability electronic_device retired unemployed housing_conditions dsf_home mobile apartment wap liheap gov_assistance_year 


** Disconnect
ologit disconnect_cat $demographics $dwelling $govassist i.state i.month [pweight=weight], vce(cl state) or
estimates store disconnect
** Table S7, Column 3
outreg2 using table_s7, cti(Disconnected from service) append bdec(2) 10pct drop(i.state) addstat(Pseudo R-squared, `e(r2_p)') eform label

** Predicted probabilities
ologit disconnect_cat age i.female i.black i.hispanic i.other i.under100fpl i.between100150fpl i.educ i.hh_under5 i.hh_disability i.electronic_device i.retired i.unemployed i.housing_conditions i.dsf_home i.mobile i.apartment i.wap i.liheap i.gov_assistance_year i.state [pweight=weight], vce(cl state) or

margins, dydx(black hispanic hh_under5 electronic_device housing_conditions wap) predict(outcome(0)) 
** Table 1, Column 5
margins, dydx(black hispanic hh_under5 electronic_device housing_conditions wap) predict(outcome(1))
** Table 1, Column 6
margins, dydx(black hispanic hh_under5 electronic_device housing_conditions wap) predict(outcome(2))

** Test of parallel lines assumption
omodel logit disconnect_cat age female black hispanic other under100fpl between100150fpl educ hh_under5 hh_disability electronic_device retired unemployed housing_conditions dsf_home mobile apartment wap liheap gov_assistance_year

** Figure 4
** Coefplot of estimates from bill, notice, and disconnect
coefplot (bill, label(Unable to pay bill)) (notice, label(Received shut-off notice) msymbol(T)) (disconnect, label(Disconnected from service) msymbol(S)), keep(age female black hispanic other under100fpl between100150fpl educ hh_under5 electronic_device retired unemployed housing_conditions dsf_home mobile apartment otherhome gov_assistance_year liheap wap) xline(1) legend(symxsize(2) row(1)) msize(tiny) ylabel(,labsize(vsmall)) xlabel(,labsize(vsmall)) xlabel(-2(2)30) legend(size(vsmall)) eform saving(figure4, replace)


** Analysis of number of instances of EI (nbreg count models)
estimates clear

** Bill pay
nbreg hhbill_number $demographics $dwelling $govassist i.state [pweight=weight], vce(cl state) irr
estimate store bill_count
** Table S8, Column 1
outreg2 using table_s8, cti(Unable to pay bill) replace bdec(2) 10pct drop(i.state) addstat(Pseudo R-squared, `e(r2_p)') eform label

** Notice
nbreg notice_number $demographics $dwelling $govassist i.state [pweight=weight], vce(cl state) irr
estimate store notice_count
** Table S8, Column 2
outreg2 using table_s8, cti(Received shut-off notice) append bdec(2) 10pct drop(i.state) addstat(Pseudo R-squared, `e(r2_p)') label eform

** Disconnect
nbreg disconnect_number $demographics $dwelling $govassist i.state [pweight=weight], vce(cl state) irr
estimate store disconnect_count
** Table S9, Column 2
outreg2 using table_s8, cti(Disconnected from service) append bdec(2) 10pct drop(i.state) addstat(Pseudo R-squared, `e(r2_p)') label eform


** Analysis of successive months 
estimates clear

** Bill
logit bill_atleast2 $demographics $dwelling $govassist [pweight=weight], vce(cl state) or
estimates store bill
** Table S9, Column 1
outreg2 using table_s9, cti() replace bdec(2) 10pct drop(i.state) addstat(Pseudo R-squared, `e(r2_p)') label eform

logit notice_atleast2 $demographics $dwelling $govassist [pweight=weight], vce(cl state) or
estimates store notice
** Table S9, Column 2
outreg2 using table_s9, cti(Received shut-off notice) append bdec(2) 10pct drop(i.state) addstat(Pseudo R-squared, `e(r2_p)') label eform

logit disconnect_atleast2 $demographics $dwelling $govassist [pweight=weight], vce(cl state) or
estimates store disconnect
** Table S9, Column 3
outreg2 using table_s9, cti(Disconnected from service) append bdec(2) 10pct drop(i.state) addstat(Pseudo R-squared, `e(r2_p)') label eform








